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We study the efficiency of algorithms simulating a system evolving with Hamil- 
tonian H = Yll^i^j- W^ consider high order splitting methods that play a key 
role in quantum Hamiltonian simulation. We obtain upper bounds on the number 
of exponentials required to approximate e~*^* with error e. Moreover, we derive the 
order of the splitting method that optimizes the cost of the resulting algorithm. We 
show significant speedups relative to previously known results. 
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^ . I. INTRODUCTION 



c 



While the computational cost of simulating many particle quantum systems using classi- 
j^ . cal computers grows exponentially with the number of particles, quantum computers have 

^ I the potential to carry out the simulation efficiently [lH3]. This property, pointed out by 

^- Feynman, is one of the fundamental ideas of the field of quantum computation. The simu- 

lation problem is also related to quantum walks and adiabatic optimization p-flOl- 
cn ■ A variety of quantum algorithms have been proposed to predict and simulate the behavior 

^ ! of different physical and chemical systems. Of particular interest are splitting methods that 

,__! ' simulate the unitary evolution e~*^*, where H is the system Hamiltonian, by a product of 

^f) . unitary operators of the form e"*"^'*', for some ti, I = 1, . . . ,N, where Ai G {-^i, • • • , Hm}, 

H = Y2"^=i ^j ^^"^ assuming the Hamiltonians Hj do not commute. It is further assumed 
VO ' that the Hj can be implemented efficiently. Throughout this paper we assume that the 

O , Hj are either Hermitian matrices or bounded Hermitian operators so that \\Hj\\ < oo for 

' j = 1, . . . , m, where || ■ || is an induced norm [T3]- 

As Nielsen and Chuang pj., p. 207] point out, the heart of quantum simulation is in the 
Lie- Trotter formula 

From this one obtains the second order approximation 

^-i(m+H2)At ^ ^-iHiAt^-iH-iAt ^0(\At\'^). 

A third order approximation is given by the Strang splitting 

^-i{Hi+H2)At ^ ^-iHi At/2 ^-iHiAt^-iHi At/2 _j_ Qn^f\3\ 

Suzuki 12, |l3j uses recursive modifications of this approximation to derive methods of order 
2k + l,ioT k = 1,2^ ••• 

A recent paper [3] shows that Suzuki's high order splitting methods can be used to derive 
bounds for the number A^ of exponentials, assuming the Hj are local Hamiltonians. These 
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bounds are expressed in terms of the evolution time t, the norm \\H\\ of the Hamiltonian H, 
the order of the sphtting method 2k + 1, the number of Hamiltonians m, and the error e in 
the approximation of e~*'^*. In this paper we will show how these bounds can be significantly 
improved. 

Consider the Hamiltonians indexed with respect to the magnitude of their norms ||-ffi|| > 
||-f^2|| > ■ ■ ■ > ll-^m||- Then the number of necessary exponentials N generally depends on 
Hi, but it must also depend explicitly on H2 since only one exponential should suffice for 
the simulation if 11-^211 -^ 0. This observation is particularly important for the simulation of 
systems in physics and chemistry. To see this, suppose m = 2 and that Hi is a discretization 
of the negative Laplacian —A, while H2 is a discretization of a uniformly bounded potential. 
Then e~*^i*i and e~*'^2<2 (^g^^ \^q implemented efficiently for any ^1,^2? and 11-^211 ^ ll-f^ill- We 
will see that, not only in this case but in general, the number of exponentials is proportional 
to both \\Hi\\ and ||i!^2||, i-e., the Hamiltonian of the second largest norm plays an important 
role. 

Let e be sufficiently small. The previously known bound for the number of exponentials, 
according to [^], is 

A^ < ATprev := m5^\m\\H\\ty+^^e-'^^^''\ (1) 

This bound does not properly reflect the dependence on H2. 

Performing a more detailed analysis of the approximation error by high order splitting 
formulas, it is possible to improve the bounds for A^ substantially. The new estimates lead 
to optimal splitting methods of significantly lower order which greatly reduces the cost of 
the algorithms. 

We now summarize our results. Recall that the Hj can be implemented efficiently but 
do not commute and Hi^iH > ||-f^2|| > ■ ■ ■ ||-f^m||- We show the following: 

1. A new bound for the number of exponentials A^, given by 



2. A speedup factor of 



N^,^ ^^^4e||iJ2r'/'' 



N ~ S'^ \ II /7i II 

-i * prev <J \ 1 1 -^ 1 1 1 

3. We show that the optimal A;*^^ that minimizes A^new is 



^new := round W-logas/g 



1, Aemt\\H2\ 

2 



On the other hand, from |^] the bound for A^prev is minimized for 

k: 




4. For k*„, the value of A^npw satisfies 



new 



2,/ii,f 1,.4-™*II^2|I 



iV^ew < 3 (2m - 1) met \\Hi\\ e^V 2 3 



For k^j.^^ the value of A^prev is 



N;^^^ = 2m'\\H\\t ■ eVi-5i-(-ll^ll*A). 



Hence 



, 1 25 4emt||H2|l / ni||Hi||t 

A^* ~ 3 ^ 

prev " 

II. SPLITTING METHODS FOR SIMULATING THE SUM OF TWO 

HAMILTONIANS 

We begin this section by discussing the simulation of 

where if i, if 2 are given Haniiltonians. Restricting the analysis to m = 2 will allow us to 
illustrate the main idea in our approach while avoiding the rather complicated notation 
needed in the general case, for m > 2. The simulation of the Schrodinger equation of a p- 
particle system, where Hi is obtained from the Laplacian operator and H2 is the potential, 
requires one to consider an evolution operator that has the form above; see [3(]. 

In the next section we deal with the more general simulation problem involving a sum of 
m Hamiltonians, Hi, ... , Hm, as Berry et al. [3] did, and we will show how to improve their 
complexity results. 

Suzuki proposed methods for decomposing exponential operators in a number of papers 
[l2 . [l3l |. For sufficiently small At, starting from the formula 

^2(^1,^2, At) = e-^^^^*/2g-^H.Aig-^/fiAt/2^ 

and proceeding recursively, Suzuki defines 

S2k{Hi,H,,At) = [S2k-2{Hi,H2,PkAt)]^S,k^2iHi,H2,{l ~ Apk)At)[S2k-2iHi,H2,PkAt)]^, 

for A; = 2, 3, ■ ■ ■ , where p^ = (4 — 4-'^/*^^^^-'^^)^-'^, and then proves that 

\\e-'iH,+H,)At _ S2k{Hi,H2,At)\\ = 0(|Atp'^+^). (2) 

Suzuki was particularly interested in the order of his method, which is 2A; + 1, and did not 
address the size of the implied asymptotic factors in the big-0 notation. However, these 
factors depend on the norms of Hi and H2 and can be very large, when Hi and ii2 do 
not commute. For instance, when Hi is obtained from the discretization of the Laplacian 
operator with mesh size h, \\Hi\\ grows as h~^. Since h = e, we get ||iii|| = 0{\). Hence, 
for fine discretizations ||ifi|| is huge, and severely affects the error bound above. 
Suppose llifill > ||ii'2||- Since 

^-i(Hi+H2)t ^ ^-iiHi+n2)\\Hi\\t 

where Tij = iij/||iii||, for j = 1,2, we can consider the simulation problem for "Hi + "Ha 
with an evolution time r = ||ifi||t. 

Unwinding the recurrence in Suzuki's construction yields 

K K 

S2k{'Hi, ■H2, At) = H S2CH1, 7^2, Z,At) = II [e-^n^z,At/2^-in2Z,At^-^n^z,At/2^ ^ ^3) 



where K = b^ ^ and each z^ is defined according to the recursive scheme, £ = 1, . . . , i^. In 
particular, zi = zk = Y\r=2Pr^ and for the intermediate values of d the z^ is a product of 
k — 1 factors and has the form Z(, = Ylreio P^ Ylreh ^^ ~ ^P^)^ where the products are over the 
index sets Jq, h defined by traversing the corresponding to i path of the recursion tree. 

Let qr = inax{pr, ^Pr — 1}, ^ > 2. Then {q,.} is a decreasing sequence of positive numbers 
and from [l4 . p. 18] we have that 

3 -A- 4A; 

Thus 

Ak 
|^,|<— forain=l,...,ir. (4) 

Equation ([3]) can be expressed in the more compact form which we use to simplify the 
notation. Namely, 

52fc('Hi, "Ha, At) = e-*^i'o^*e-^^2"^^*e-^^^'i^* ■ ■ ■ g-^^^^^Aig-i^i.^At^ ^5^ 

where Sq = -2i/2, Sj = {zj + Zj+i)/2, j = 1, . . . ,K — 1, and sk = z^/^. Observe that 
l^j=0^j ~ -*-' 2^j=i ^j ~ ^■ 

K 



We need to bound 0"^ = 'Y1,j=q \^j\ + ^j=i l^il from above. From (jl]) we have 



and also 



/ J I Jl — ok ' 

.=1 "^ 



^ 4A;5'=-i 



j=0 

Thus 

0-fc < -A; ( - j =: Cfc for fc > 1. (6) 

(The above trivially holds for k = 1.) 

Expanding each exponential in ([Sj) we obtain 

S2kini,n2, At) =(/ + -Hisoi-iAt) + ^nlsli-zAtf + ■■■ + ln1s',{-zAt)' + ■■■) 
(/ + ■H2Zi{-tAt) + ^nlzii-iAty + ■ ■ ■ + ^n'^zli-iAt)'' + ■■■) 

(/ + nis,i~iAt) + ^nlsli-iAtf + ■■■ + l-H^st(-^At)'= + ■ ■ ■ ) ^^^ 

(/ + n2ZK{-iAt) + hilzl^-iAtf + ■■■ + ^•H2'4(-^At)'= + ■ ■ ■ ) 
2 k\ 

(/ + HiSK{-iAt) + ]^Hlsl{~iAtf + ■■■ + i-Hj4(-^At)^ + ■■■)• 



After carrying out the multiplications we see that 5*2^ is a sum of terms that has the form 

aolail ■ ■ ■ axU^i^- ■ ■ ■ Pk^- 

where the ao,ai,- ■ ■ , ax and the /3i, ■ ■ ■ , /3k are obtained by multiplying the denominators 
in the expansion of the exponentials. 

The terms that do not contain 7^2 are those for which Pi = (32 = ■ ■ ■ = (3k = 0, and their 
sum is 

^^-^ ao\ai\---aK\ 






(9) 

■1 qa- 

X K 

= Y[e-'^''^^' = exp(-z^His,At) = exp(-iHiAt). 
i=o j=o 

On the other hand, consider 

^.,^m+H,)At ^j^ (_^(^^ ^ ^^)^^) ^ . . . ^ :^(-'C^i + ■^2)At)'= + ■ • • . (10) 

The terms that do not contain 7^2 sum to 

00 _ 
Y,^HiH^t)'' = e-'^^^'. (11) 

fc=0 

Let us now consider the bound in ([2]). Clearly the terms that do not contain 7^2 cancel 
out. Therefore, the error is proportional to ||'H2|||Atp^"'"^, i.e. it depends on the ratio 
||-f^2||/||-f^i|| of the norms of the original Hamiltonians. This fact will be used to improve the 
error and complexity results of Berry et al. [3] 

Lemma 1. For A; e N, Cfe|At| < k + I (see, Eq. \B^ and \\'H2\\ < ||'Hi|| = 1 we have 

II exp(-z(Hi + ?^2)At)- 52^(^1,^2, At) II < JiM^(cfc|At|)2'=+i. (12) 

[Zk + 1)! 

Proof. For notational convenience we use 5'2fe(At) to denote S2k{(Hi,'H2,At). Consider 

00 
exp{-t{n,+n2)At)-S2kiAt)= J2 [Ri{^t)-Ti{At)], (13) 

l=2k+l 

where Ri{At) is the sum of all terms in exp(— 'j('Hi+'H2)At) corresponding to At' and Ti{At) 
is the sum of all terms in 6*2^ (At) corresponding to At . Moreover, we know that the terms 
with only "Hi cancel out. Hence, we can ignore the terms in T;(At) and Ri{At) that contain 
only Til (and not 7^2) as a factor. It follows that 

Ri{At) = 1(7^1 + 'H2yHAty - ^H[{-iAt)K (14) 



Then 

\\RiiAt)\\<^2'\\n2mt\\ (15) 

since there are 2^ — 1 terms, and they are bounded by jf||'H2|||^^|^- 
Now consider the terms in Ti{At). From f l7|8|) 






(16) 
where the condition Yli=i A 7^ hold because there are no terms containing Tii alone. Since 
the norm of W^^H^^'UT ■ ■ ■ H^'Hi" is at most ||-H2||, we have 

Note that we relaxed the condition X]i=i A 7^ since it does not affect the inequality. 

To calculate the sum V ° .^ , — ^ .„\ „-^, , where V- „ ^j + t^,- t Pi = ^^ we first consider 
the following equation 

exp(|soAt|) exp(|2;iAt|) exp(|siAt|) ■ ■ ■exp(|2;i^At|) exp(|s/<At|) 

.ao=0 ^ / \/3i=0 '^ / \ai=0 






/3x=0 / \aK=0 



aQ\ai\---aK\Pi\--- Pk^ 



p=0E°.+E/3i=P 



Hence Ky-^ , v-o , " ,\ — ^ .} , J^, is the coefficient of At in the equation above. Sim- 

^^^E"j+EPj=' ao!ai!---ax!/3i!---/3K! ' ' ^ 

ilarly, 



exp(|soAt|) exp(|ziAt|) exp(|siAt|) • ■ ■Qx:^{\zxAt\) exp(|sif At 

K K 

= exp((^ |s,| + J2 kil)|At|) = exp(afe|At|) 



K K 

'UiHIAtn =exD('CTjAtn 

(19) 



i=0 1=1 

oo 



p=0 •'^ 



Recall that the bound for cr^ given in Eq. ([6]). Thus the coefficient of |At|' is bounded from 
above by ^c\. Therefore, we have 

||TKAt)||<^||-H2|||At|^ (20) 



We combine Eq. (^, ([20]), to obtain 

oo 

\\exp{{ni+n2)At)-S2k{At)\\< J2 II^KAi)-T,(At)|| 



l=2k+l 



< J2 \\Ri{At)\\ + \m{At)\\ 



l=2k+l 



oo ; 



<2 y: "f\m\At 

l=2k+l 

2 



(21) 



< 



< 



{2k + l)\ 

4 
(2A; + 1)! 



1^211 |CfcAt| 



2k+l 



2A; + 2 



-1 



\n2\\\ckAt\''+\, 



where the last two inequahties follow from the assumption Cfc|At| < k + 1. and an estimate 
of the tail of the Poisson distribution; see, e.g., [iSl, Thm 1]. D 

Theorem 1. Let 1 > e > be such that 8et\\H2\\ > e. The number N of exponentials for 
the simulation of e^^^^^^^'^^^ with accuracy e is bounded as follows 



N <3 5^-^ 



\Hi\\t 



8et\\H2\\V^^^^'^ 8e /5' '"'^ 



3 V3 



for any A; G N, where \\H2\\ < ||-f^i||- 

Proof. Let M = lAtl^"*^. Then using Lemma 1 and Tij = Hj/\\Hi\\, j = 1,2, we obtain 



|,-.(^.+...)*_^M||H.||*^^^^^^^^/^)|| ^^ll^^ll^ 4_^ll^^ll (^) 



2k+l 



mH2 



{2k + l)\ 

„2fc+l -I 



{2k + l)\M^f'' 
Recall that Ck is defined in ([6]) and is used in Lemma 1. For accuracy e we obtain 

-\e{2k + l)\J ■ 

We use Stirling's formula jl6|, p. 257] for the factorial function 

{2k + 1)! = V2^i2k + i)(2fc+i)+i/2e-{2fc+i)+e/(i2(2fc+i))^ < ^ < 1, 
which yields 



[(2A; + l)!]-i/(2fc)<gi+i/(2fc)/(2A; + i). 



(22) 



It is easy to check that 
Thus it suffices to take 



l/(2fc) ^ 9l+l/{2fc) 



M> 



k 

8et\\H2 



l/{2fe) 



2e Cfc 
2A; + 1' 



So we define M to be lower bound of tlie expression above, i.e., 



M:-- 



Set\\H2\ 



l/(2fc) 



2ecfc 
2k + l' 



It is easy to clieck tliat 



2e 



2k + I 



{k + l) >e. 



wliicli along with the condition 8et||if2|| > ^ yields M{k + 1) > Ck- This shows the assump- 
tions of Lemma 1 are satisfied with this value of M. 

From the recurrence relation the number of required exponentials to implement 5*2^ in 
one subinterval is no more than 3 ■ 5^~^. We need to consider two cases concerning M||ifi||t. 
If M||ifi||t > 1, then the number of subintervals is [M||iiri||t], i.e., we partition the entire 
time interval into an integer number of subintervals, each of length at most M~^. The total 
number of required exponentials is bounded by 3 • 5*^~^[M||iJi||t]. Substituting the values 
of M and Ck we obtain the bound for A^. In particular. 



N <?>-h 



fc-i 



\Hi\\t 



SetlliJsl 



l/(2fc) 



8e 

y 



k-i 



(23) 



If M||i^i||t < 1, then Lemma 1 can be used with At = ||-f/^i||t, since ||ii/^i||t < M ^ and 
we have already seen that M is such that the assumptions of Lemma 1 are satisfied. Thus 



|e-*(^^+^^)*- ^2.(^1,^2, ||i/i||t)|| < 



{2k + I] 



1^211 (C,.||f^l||t) 



2k+l 



At\\H2 



c 



2k+l 



{2k + I) 



HiW < At\\H2 



c 



2k+l 



i2k + r 



AMY'' < e. 



where the last inequality holds by definition of M. In this case the total number of expo- 
nentials is simply 

A^ < 3 ■ 5''-\ (24) 

Combining ( 123|) and (!24l) we obtain 



N <3-5''-^ 



\Hi\\t 



8et||i72| 



ym g^ .5. k-1 



This completes the proof. 



D 



Remark 1. Lemma 1 and Theorem, 1 indicate that when \\H2\\t <^ e then the number of 
exponentials N can he further improved. In this case it can he shown that high order splitting 
methods may lose their advantage. We do not pursue this direction in this paper since we 



assume that the Hj, j 



, m, are fixed and study N as e -^ 0. 



III. SPLITTING METHODS FOR SIMULATING THE SUM OF MANY 

HAMILTONIANS 

In this section we deal with the simulation of 



-ij:T=iHjt 



where Hj, j = l,...,m, are given non-commuting Hamiltonians. The analysis and the 
conclusions are similar to those of the previous section where m = 2, but the proofs are 
much more complicated and certainly tedious. This is the problem that Berry et al. [3] 
considered. 

We use Suzuki's recursive construction once more [l3|. In particular, for 



j=l j=m 

and 

S2kiHi, . . .,Hm,At) = [S2k-2{PkAt)]^S2k-2{{l " 4pfc) At) [^2fc-2(Pfe At)]^, A; = 2, 3, ... , 

where for notational convenience we have used S2k~2{'At) to denote S2k-2{Hi, ■ ■ ■ , Hm, At), 
and pfc = (4 - 4^/^^''"^))"^ we have that 

ll^-.E- iH,A* _ g^^^jj^^ jj^^ ^^)|| ^ o(|Ai|2fc+i)_ (25) 

Assuming again that \\Hi\\ > \\H2\\ > ■■■ > ||-f^m|| we normalize the Hamiltonians by 
setting Tij = Hj/\\Hi\\, j = 1, . . . ,m, and consider the equivalent simulation problem 

where r = ||-f/^i||t. Proceeding in a way similar to that for m = 2 of the previous section we 
derive the following lemma, whose proof can be found in the Appendix. 

Lemma 2. For k eN, 4|At| <k + l,dk = m{A/3)k{5/3)''-^ and \\Hm\\ < ■■ < ||^2|| < 
llT/ill = 1 we have 



4||-H2 
■Hj/\t) - ^2k(:Hi, ..., -Hm, i\t) II < - 

i=i 
From Lemma 2, we have the following theorem. 



expH5]?^,At)-^2.(^i,...,^™,At)|| < ^^{dk\At\f^+\ (26) 



Theorem 2. Let 1 > e > he such that 4met||iJ2|| > ^- The number N of exponentials for 
the simulation of e^^^^^^ \-Hm)t yjj^ifi accuracy e is hounded by 



l/(2fe) . „ / r:\k~l 



/ 4emt||i72|i y^^ Mme /5\ 



A^< (2m- 1)5'="^ 

for any fc G N, where ||-f^m|| ^ ■ ■ " ^ ll-^2|| ^ ll-^ill- 

Proof. The proof is similar to that of Theorem 1. Let M = |At|^^. Then using Lemma 2 
and Tij = Hj/\\Hi\\, j = 1, . . . , m, we obtain 



\^-i(Hi+-+H^)t nM\\Hi\\t,^ ^ 1/A/r^|| < MWH \\f Wl-I II (^\ 

|e - ^2k [f-Li,..., Hm,i-/M)\\ S M||/ii||t _^ -|^N,II^2|| 1^1 



2fc+l 



r2k + l . 

nH2\' "■ 



(2A; + l)!M2'=' 



10 



Recall that dk is defined in Lemma 2. For accuracy e we obtain 



M> 



4^11^2 II 4 



2k+l\ l/(2fc) 



e{2k + l)\ 
We use the estimate ( 122|) . It is easy to check that 

Thus it suffices to take 

/4emt||iJ2||\^^^^^'^ 2e4 



M> 



V £ / 2fc+l 

So we define M to be the lower bound of the expression above, i.e., 

nemt\\H2\\V^^^''^ 2edk 



M ■-- 



/4emt||i72|| Y 



2fc + l' 



As in the proof of Theorem 1, it is straightforward to verify that M{k + 1) > dj^. Therefore, 
the assumptions of Lemma 2 are satisfied for this value of M. 

From the recurrence relation, we see that the number of required exponentials to imple- 
ment 5*2*; in one subinterval is no more than (2m — 1) ■ 5^^^. Again we distinguish two cases 
for M||ifi||t. We deal with the case M||i!fi||t < 1 in the same way we did in the proof of 
Theorem 1, to conclude 

N <{2m-l)- 5''-\ 

If M||i7i||t > 1, then the total number of required exponentials is 

A^< (2m-l)-5'^"^[M||ifi||t]. 
Substituting the values of M and dk we obtain 

V(2fc) A^^ /K\k-l 



N <{2m-l)-5 



fe-i 



/'Aemt\\H2\\\ ^^ ' Ame fbY 



This completes the proof. D 

The reader may wish to recall Remark 1 that applies in the case too. 

Corollary 1. If in addition to the assumptions of Theorem 2 either of the following two 
conditions holds: 

• 4met\\Hi\\ > 3 

• e is sufficiently small such that 

, 4met||i/i||\ „, 5, 4met||iJ2|| 

m — 2 In - In < 

5 / 3 e 

then the number of exponentials, N, for the simulation o/e^*'-^^^"^''^'"^* with accuracy e is 
hounded by 

, .1,, „ fAemtWH^WV'^^'^Hmef^V-^ 
N<2i2m-l),^-^\\H4t[-^) ^(3) , 

for any A; G N. 



11 

Proof. From the assumption of Theorem 2 we have 4emt\\H2\\/e > 1. Consider the argument 
of the ceihng function in the bound of Theorem 2. It is greater than or equal to 1, if 
4met||ifi|| > 3. Otherwise, we take its logarithm and multiply the resulting expression by 
k. This gives the quadratic polynomial 

2, 5 , ^, , Amet\\Hi\\ ^ , 4met||i;f2|' 



2k^ \n- + 2k In ^^ + In ■ 



3 5 e 

When e is sufficiently small and the discriminant is negative, i.e., when 

In ^^ - 2 In - In ^^ < 0, 

5 / 3 e 

the polynomial is positive for all k. Hence, that argument of the ceiling function in the 
bound of Theorem 2 is greater than 1, for all A; > 1. 

In either case, we use [x] < 2x, for x > 1, to estimate A^ from above. D 

IV. SPEEDUP 

Let us now deal with the cost for simulating the evolution e~*'-^j=i^^^*. Berry et al. Q 
show upper and lower bounds for the number of required exponentials. We concentrate on 
upper bounds and improve the estimates of [^. 

We are interested in the number of exponentials required by the splitting formula that 
approximates the evolution with accuracy e. Recall that 

iV„:^2,2»-l,5-||...|K(lH«)""'l|£(|)'" 

exponentials suffice for error e. The above estimate holds for e sufficiently small as Theorem 2 
and Corollary 1 indicate. The corresponding previously known estimate [3] is 

j_ 

, / I \ 2fc 

Np,,, = m 5^" im\\Hi\\ty+-^ 

where H = ^ -^^ Hj. 

The ratio of the two estimates is 

"".4fii#V". (27) 



A^prev ~ 3H 11^1 . 

So for large k we have an improvement in the estimate of the cost of the algorithm. On 
the other hand, if ||-f^2|| ^ ll-^^ill we have an improvement in the estimate of the cost the 
algorithm not just for large k but for all k. This is particularly significant when k is small. 
For instance, k = 1 for the Strang splitting 5*2, which is frequently used in the literature. 

Let us now consider the optimal k, i.e., the one minimizing N^ew, for a given accuracy e. 
It is obtaind from the solution of the equation 

2, 25 ^^Aemt\\H2f 



2k^ In — - In 



3 e 



12 



Since we seek a positive integer A;*^^ minimizing A^new, "we set 



1 4:emt\\H2\ 

2 



4w := max <^ round \ - logss/g Z > 1 f , 



where round(x) = \_x + 1/2J, x > 0. We can avoid using the max function in the expression 
above by considering e < mt\\H2\\. Then the number of exponentials A^new satisfies 



„ , 1 , 25 , 4emt\\H2\\ 
2\/oln-:rln- 



iV^ew < 3 (2m - 1) met \\Hi\\ e^ 2 3 



Berry et al. Q find 



1 /, m\\Hi\\t 



k;^,^ = round I - \l logs ^T^ + 1 h (28) 



which minimizes A^pre- For k*^^ the number of exponentials A^prev becomes 



'prev 

As a final comparison with A^prev we have 



mll/fillt 

AT* = 2m^H,\\t e'V^'^^'°^— . (29) 



, ;i 25 4emi||_H'2|| / mllHillt 

N* ~ 3 '^ 

prev " 

Hence, there is an important difference between the previously derived optimal k and the 
one derived in the present paper. In Q], the optimal k depends on ||-ffi||. More precisely, 
we show that the optimal k depends on ||-ff2||, the second largest norm of the Hamiltonians 
comprising H, which can be considerably smaller than ||-f^i||. 
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VI. APPENDIX 

Proof of Lemma 2. Unwinding the recurrence for S2k we see that 

K K 



S2k{T^l-, ■ ■ ■ , 'Hm, At) — J_J_5'2(('Hl, . . . , Tim, ^lAt) — J^J^ 



£=l i=l lj=l j=m 



J-Jg-.H,.,Ai/2j-Jg 



■iHjZiAt/2 



where K = 5'' ^ and each zg is defined according to the recursive scheme, i = 1, . . . ,K. For 
the details, see the part of the text that follows ([3]). The bound (j4]), namely, 

Ak 
\ze\ < -^ ioT a\\i=l,...,K, 



13 

holds independently of m, because it depends on the k — 1st levels of the recursion tree and 
not on the leaf, 5*2 ((Hi, . . . , Tim, ziAt), at which, the corresponding to i, path ends. 

In the expression of S2{{'Hi, ■ ■ ■ yT-Lm, ziAt) the sum of the magnitudes of the factors 
multiplying the Hamiltonials in the exponents is m\ze\ ■ |At|, for all i = 1,...,K. Thus 
in the expression of 5*2^ above, the sum of the magnitudes of all factors multiplying the 
Hamiltonians in the exponents is 

A 4k 

^(ml^^l • |At|) < S^-^m— |At|. 

Define 

4 := mh ("0 A; > 1. (30) 

Equivalently, one can view the expression for 5*2^ above as a product of exponentials of 

the form e^^*"^-"^*, where YlnLi fj,n = 1; J = I5 ■ ' ' 5 ^^5 ^ind Nj is the number of occurrences 
of l-Lj in S2k- Recall that for m = 2 we used s„ to denote ri „ and Zn to denote r2,„. With 
this notation and using fl30|) we have 



J^ki,n| < 4- (31) 



Jl" 



(Recall the derivation of ([6]).) 

Expanding the factors of 6*2^ in a power series individually, and then carrying out the 
multiplications amongst them, we conclude that S2k is given by an infinite sum whose terms 
have the form 

n-^Hj-Hr,,.At]^-. (32) 

The factors of these products are specified by the Hamiltonians Hj and the order of their 
occurrences after unwinding the recurrence for 6*2^, where j = 1, . . . , m and 7j,„ = 0, 1, 2, ... , 
for all n = 1, . . . , Nj. 

Consider the terms that contain only T-ii and, therefore, have 7j_„ = 0, for ra = 1, . . . , Nj 
and j = 2, . . . ,m. The sum of these terms is 



1 °° 1 

7i,i=---=7i,iVi=0(l,n) '^^'"' 



7i,n=0 for jVI 0» '■^'"' 7i_i=...=7i_jVi=0(l,n) '^ 



n=l 7i,n '^^''''' n=l 

= g-*E„''-i,nHiAi ^ g-J-HiAt 



(33) 



On the other hand. 



-zJ]?^,Atj+--- + -f-z5^H,Atj 



(34) 
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and the terms that contain only "Hi have sum 



At 



fc=0 



(35) 



Let us now consider the error bound in (!25|) . The sum of the terms with only T^i in 
S2k+i and G^p{YlT=i HjAt) is the same and cancels out when we subtract one from the 
other. Moreover, in exp(— -i ^"1^ "HjAt) — S'2fc(At) we know that the terms of order up to 
2k also cancel out, see Eq. (!25|) . From this we conclude that the error is proportional to 

11-1/ |||Aj-|2fc+l 

Consider 

OO 

eM-^i'Hi + --- + 'H„^)At)-S2k{'Hi,...,'Hm,At)= J^ [M^t) -Ti{At)], (36) 

l=2k+l 

where Ri{At) is the sum of all terms in exp{~i{'Hi + ■ ■ ■ + 'Hrn)At) corresponding to At' 
and Ti{At) is the sum of all terms in 5*2^; corresponding to At . We can ignore the terms in 
Ti{At) and Ri{At) that contain only Hi (and not I-L2) as a factor. 
Then 



\R,{^m 



i(|:«,A<)-i«; 



ttH At' 



<— ||?/2|||At| 



(37) 



because there are m'- — 1 terms in Ri and each norm is at most i||'H2lllAt|'. 



From (1321) we have 



IT II "-2 1 



TiiAt)= Yl 









Yl nf'-At^ 



(38) 



where ^„7i,n 7^ I, i-e., there is no terms containing only Hi. So, || Hon) '^/'"ll — ll'^2 
and 

||TKAt)||< Yl fT ' , ll^2|||Atr. 



(39) 



E7j,n 

To calculate the coefficients of the sum, we consider 



00 _ 
nexp(|r,,„At|) = J] Y :^7l^^-.«^^r" 



(i,n) 



(i.™)7j,n=0 

00 



7i,. 



nir- I'^J'" 

z^ n -Y- ' ' ' 

'=0 E7,>=« ^^^> '^'"' 

Hence the coefficient of |At|' in ( 1391) is equal to that in (HUj) . Also 



(40) 



]^exp(|rj,„At|) = exp(^ |rj'„At|). 



J,n 



],n 



From (13T!) we obtain 



|TKAt)|| = ^||?^2|||At|^ 



(41) 



(42) 
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Therefore, 



expiJ2'^Ai) - S2kiAt)\\ < J2 \\Rii^t)\\ + \\TiiAt)\\ 

j=l l=2k+l 



oo 



<2 J2 7fll^2|||At|^ 



l\ 

l=2k+l 



°° 1 

2\m\ y: hdkAt\^ (^^) 



n 

l=2k+l 



<7^J-^\\n2\\\dkAt\^'^'(i--'^'^^'^^" 



(2A; + 1)!" '"' " ' V 2A; + 2 

where the last two inequahties follow from the assumption dk\At\ < k + 1 and an estimate 
of the tail of the Poisson distribution; see, e.g., [11, Thm 1]. D 
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